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ABSTRACT 


A rarefied gas with hard- sphere molecules enclosed between parallel walls is ana- 
lyzed for the case of couette flow and heat transfer by a model sampling procedure. In 
one case, both walls are stationary and at different temperatures. In the other case, the 
upper wall is moving in a direction parallel to the lower wall and both walls are at the 
same temperature. The wall accommodation factor is taken as unity. The target mole- 
cule velocity distribution is presumed to be the sum of two different half-Maxwellians. 
The parameters describing the half-Maxwellians are assumed constant over a given zone. 
By scoring the properties of sample molecules as they pass scoring positions, the mac- 
roscopic quantities, of interest are obtained. The problem is iterated until, these values 
agree with the assumed distribution values. 



RAREFIED-GAS COUETTE FLOW AND HEAT TRANSFER 
BETWEEN PARALLEL PLATES BY MODEL SAMPLING 
by Morris Perlmutter 
Lewis Research Center 

SUMMARY 

A rarefied gas with hard- sphere molecules enclosed between parallel walls is ana- 
lyzed for the case of couette flow and heat transfer by a model sampling procedure. In 
one case, both walls are stationary and at different temperatures. In the other case, the 
upper wall is moving in a direction parallel to the lower wall and both walls are at the 
same temperature. The wall accommodation factor is taken as unity. The wall tempera- 
ture ratio and wall velocity are taken large so that linearized solutions are not applicable. 
The channel between the walls has been divided into zones with boundary surfaces paral- 
lel to the walls. The target molecule velocity distribution in each zone is presumed to be 
the sum of two different half-Maxwellians. The parameters describing the half- 
Maxwellians are assumed constant over a given zone. By scoring the properties of sam- 
ple molecules as they pass scoring positions, the macroscopic quantities of interest, 
such as temperature, density, shear, and heat transfer are obtained. Also, the average 
thermal velocity and number density in each zone are found and compared with the values 
assumed for the target molecule distribution in that zone. The problem is iterated until 
these values agree. 


INTRODUCTION 

In rarefied gas transport problems the usual macroscopic approach using Newton’s 
law of friction or Fourier’s equation for conduction is no longer applicable because these 
equations are derived assuming small changes of fluid properties over the mean free path 
lengths of the molecules. Therefore, a microscopic approach based on the kinetic theory 
of gases must be used. 

The present analysis treats a system consisting of a rarefied gas enclosed between 
two walls. Two cases are treated: in one, both walls are stationary and at different 


temperatures; while in the other, both walls are at the same temperature but one wall is 
moving. In the present model, the molecules are assumed to have hard- sphere colli- 
sions and a wall accommodation coefficient of 1. The problem is analyzed by a Monte 
Carlo technique, that is, a model sampling procedure. The Monte Carlo method has 
been extensively used in the field of neutron transport theory. Reference 1 discusses the 
Monte Carlo method and its various techniques and applications. 

Several analytical treatments of the eouette flow problem have been carried out 
(refs. 2 to 5). In general, these consisted of setting up systems of moment equations 
which are solved to obtain the space dependence of the parameters occurring in the veloc- 
ity distribution functions. Gross and Ziering used the linearized form of the Boltzmann 
equation to consider the case where the wall velocities were relatively low and both walls 
were at the same temperature (ref. 2), and the case where both walls were stationary but 
at different temperatures (ref. 3). The local molecular velocity distribution was taken 
as two-sided, each side being approximated by a different polynomial, and then half- 
range moments were formed from the linearized form of the Boltzmann equations. Liu 
and Lees (ref. 4) assumed a two-sided Maxwellian velocity distribution and Maxwellian 
molecule collisions, that is, the repulsive force between the molecules varied inversely 
as the fifth power of the distance between molecule centers. The moment equations led 
to the determination of macroscopic quantities of interest such as temperature and heat 
transfer. This same approach was used by Lavin and Haviland (ref. 5) for the case of 
stationary walls at different temperatures assuming hard- sphere interaction potentials 
for the molecular collisions. 

In the method of solution based on the Monte Carlo method, sample molecule histo- 
ries are generated by randomly picking choices at points of decision in the history of the 
sample molecule from the appropriate probability distributions. The sample molecules 
are followed through the mathematical model of the system, and by averaging molecular 
properties at various positions the macroscopic quantities describing the system can be 
obtained. This method reduces the complexity of the analysis since the relations in- 
volved the setting up a single sample molecule history are relatively simple; thus, sim- 
plifying approximations made in the analytical procedures are avoided. However, the 
Monte Carlo method requires extensive calculations on a high-speed computer to obtain 
the many sample histories needed to obtain necessary accuracy. 

The Monte Carlo method was used by Perlmutter (ref. 6) to treat a collisionless 
plasma flowing through a channel with an imposed magnetic field; by Haviland (ref. 7) to 
treat a rarefied gas enclosed between stationary parallel plates at different tempera- 
tures; and by Perlmutter (refs. 8 and 9) to determine the heat transfer through a gas be- 
tween stationary walls at different temperatures. In the analysis by Haviland (ref. 7), 
the amounts of time the sample molecule spent in the zone into which the channel had 
been divided were noted, and from these values a numerical table could be made up for 
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the molecular velocity distribution in each zone. This distribution constituted the target 
molecule distribution in the next iteration. After three or four iterations, although full 
convergence was not obtained, no evidence of the results being nonconvergent was found. 
The moments of the distribution for each zone corresponding to the macroscopic quantity 
of interest was then found by integration of the last obtained velocity distribution. In the 
analysis by Perlmutter (refs. 8 and 9) a full Maxwellian was assumed for the form of the 
target molecule velocity distribution with assumed initial values of the parameters under 
investigation. Then by summing properties of the sample molecules as they cross the 
zone boundaries, the macroscopic local conditions could be found. These new param- 
eters of the target molecule distribution served as the basis for the next iteration of the 
target molecule distribution as in Haviland’s case. Three or four iterations were needed 
to obtain the required agreement. 

The present analysis also follows a sample molecule traveling through the model. 
However, in the present case, the target molecule velocity distribution in each zone 
within the channel is assumed to be two half-Maxwellian distributions of the form as- 
sumed by Liu and Lees in reference 4. By averaging properties of the sample molecules 
at scoring position across the channel, the parameters of the distribution as well as the 
macroscopic flow properties such as temperature, velocity profile, heat transfer, and 
shear stress can be found. 
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SYMBOLS 


average thermal velocity (2R T)^^ 

© 

specific heat at constant pressure 
channel width 


error function (2/Vir) 



du 


probability density function 

quantity obtained by Monte Carlo procedure (appendix F) 

Boltzmann’s constant 

Knudsen number, X^. /D = m/\/2 p^Dttct^ 

number of values of g obtained (appendix F) 

dimensionless velocity, u/C w q 

mass of molecule 
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N 

number of sample molecules emitted from surface 

wO 


proportional to flux of molecules leaving surface 

wO 

P 

pressure, P-. /3 


p « 

shear stress, -p(VjVj) 


Pr 

Prandtl number, pc /K 

Jr 


P 

zone or scoring-position number 


Pf 

last scoring-position situated at surface wl 


Q 

property of sample molecule 


<Q) 

average quantity J' Qf d^V 


<Q) + 

average quantity J' Qf + d^V 


q 

heat-transfer rate 


q T 

total energy transfer rate 


R 

random number between 0 and 1 


R g 

gas constant, K/m 


s 

2 

mutual collision cross section ira 


s + ,s_ 

number of sample molecules through scoring position in 
Xg- direction 

T 

absolute temperature 


t 

dimensionless temperatures, T/T w q 


u i 

J.!. 

mean velocity, i component 


U i + ’ U i- 

i in mean velocity of molecules in positive, negative 

x 2 - ( 

u wO> u wl 

velocity of walls 0 and 1 


V i 

1.1. 

molecular velocity, i 11 component 


V i 

thermal velocity, Vj - u. 



transformed i Ln component of thermal velocity 


v* 

i 

velocity after collision 


V R 

relative velocity 


X i 

til 

i tn coordinate 


X. 

1 

1.1. 

dimensionless i coordinate, x./D 


r 

defined by eqs. (C6) and (C7) 
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e defined by eqs. (C6) and (C7) 

77 viscosity coefficient 

k thermal conductivity 

0 collision rate of sample molecule with target molecules 

9 angle 

X path length to collision 

Ajyj mean free path length through Maxwellian gas with density p 
locity C 

A„ mean free path length 

O 

p dimensionless velocity, V’/C 

/j. r relative dimensionless velocity, /C 

Pr a defined by eq. (BIO) 

p mass density 

a diameter of hard sphere molecule 

cp angle 

Subscripts: 

A averaged across channel 

M Maxwellian 

M=0 both walls stationary 

p zone or scoring position number 

Pj last zone or scoring position number, situated at surface wl 

R relative velocity 

r radial coordinate in cylindrical system 

s sample molecule 

t target molecule 

wO wall 0 

wl wall 1 

a, /3 subscripts used in appendix B 

0; 1 evaluated in gas adjacent to wall 0 or 1 

1, 2, 3 coordinate directions 


and thermal ve- 
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+, - positive or negative Xg-direction 

°° collisionless solution 


ANALYSIS 

The problem under consideration is the determination of the properties of a rarefied 
gas contained between two parallel walls, each infinite in area. As can be seen from fig- 
ure 1, which illustrates the geometry involved, the normal to the walls is in the x 2 ~ 
direction. One wall is located at Xg = 0, the other at Xg = 1. The wall at Xg = 0 is 
stationary and at temperature T w0 ; the other wall at Xg = 1 is moving with a velocity 
u^j in the x^-direction and its temperature is T w ^. 

The channel width is divided into zones, with boundary surfaces parallel to the walls. 
The zones act to break up the distance between walls into subintervals. Sample molecule 
histories are generated by starting out from the wall a molecule picked from a Maxwel- 
lian velocity distribution corresponding to the wall temperature and following the mole- 
cule through the various zones where it can collide with target molecules distributed 
throughout the space between the walls. In each zone, the velocity distribution of the 
target molecules is taken as uniform over the zone and consists of two half-Maxwellians, 
each half corresponding to a different temperature, density, and velocity. Thus the 
velocity distribution of target molecules per unit velocity space is given by 


^ f t 



for V 2 > 0 


> 


for V 2 < 0 


( 1 ) 


where the subscript + refers to molecules moving in the positive V 2 direction, that is, 
toward the wall wl at x 2 = 1 and the subscript - refers to molecules moving toward 
wO at x 2 = 0. It is assumed that the molecules reflected from the walls are perfectly 
accommodated, that is, the reflected molecules are assumed to have been emitted from a 
Maxwellian distribution corresponding to the wall temperature. The velocity of the sam- 
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pie molecule leaving the wall (computer flow chart position 1, see fig. 2) as derived in 
appendix A is given by 


v 2w = - c ; in r v . 


, 1/2 


(2a) 


/ 2 W 2 

v 1w = r c w ln R v r ) cos 2nR e + U W 


(2b) 


\l/2 


V 3w = \- C w ln R vJ sin 2?rR 
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(2c) 


where R is a random number between 0 and 1 chosen for each sample molecule leaving 
the surface. 


Path Length to Collision 


After the sample molecule leaves the wall it is necessary to calculate its path length 
to first collision (flow chart position 2) to find if it has undergone a collision before pass- 
ing through the first zone. The probability density function that a sample molecule will 
collide in the incremental path length A to A + dA is given (refs. 9 and 10) as 



( 3 ) 


where X is the mean free path to collision of the sample molecule moving at velocity 
s 

V in that zone. A path length to first collision for the sample molecule can be chosen 
s 

from this distribution by the same procedure used previously, that is, employing the 
random number with range between 0 and 1 and relating it to mean free path X by 
the equation 


X = -A s ln R^ (4) 

In order to pick from this distribution, the mean free path A for the sample molecule 
in the zone must be known. If the undirected velocities of the target molecules were in a 
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Figure 2. - Monte Carlo flow chart (Blocks designated positions 1 to 5 are discussed in text ! 




































Maxwellian distribution with an average thermal velocity C and density p, the mean 
free path for the sample molecule X M as shown in appendix B (eq. (12)), would be given 
by 



where Kn is the usual definition of Knudsen number for a gas m a Maxwellian distribu- 
tion with hard- sphere molecule collisions. 

However, the target molecule velocity distribution has been assumed to be two half- 
Maxwellians (eq. (1)). The direct calculation for a collision for a sample molecule with 
a target molecule picked from this distribution would be very complex. One method of 
treating the collision that introduces no additional assumptions and allows the calcula- 
tions to be greatly simplified, although at the expense of additional calculations, is to re- 
sort to the following artifice. 

The assumption is made that the distribution of target molecules consists of two full 
Maxwellian distributions: 


eh = 2p + f M+ + 2 p -*m- <6) 

where f M+ and f M _ are Maxwellian distributions with densities 2 p + and 2 p_ and 
thermal velocities of C + and C_, respectively. (It is realized that the molecules cor- 
responding to f M+ with V 2 < 0 and the molecules corresponding to f M _ with Vg > 0 
are in actuality not present in the present model. ) This assumption permits the calcula- 
tion of the mean free path using the Maxwellian distribution relations (eq. (B13) of appen- 
dix B) . Thus, the mean free path for a sample molecule moving through the two full 
Maxwellian gases is given by 


X s X M+ X M- 

where X M+ and X M are determined by equation (5), appropriately modified by + and 
- subscripts, respectively. 

The mean free path X_ calculated in this manner will be shorter than for the case 

s 

for the assumed two half- Maxwellian target distributions (eq. (1)), since collisions with 
the nonexistent target molecules are now included. These collisions with nonexistent 


10 



target molecules are made to have no effect on the sample molecule, by taking the path 
and velocity of the sample molecule as unchanged after a collision with a nonexistent tar- 
get molecule. In essence, then, the nonexistent target molecules do not interact with the 
sample molecule. 

After a path length to collision is picked for the sample molecule from equation (4), 
it is compared with the distance the sample molecule must travel to cross the zone (flow 
chart position 3). If X is greater than this distance, the sample molecule is started at 
the beginning of the next zone with its velocity components unchanged and the procedure 
as before. If it is not, there is a collision in the zone and new velocity components must 
be calculated for the sample molecule after collision. 


Sample Molecule Velocity After Collision 


To find the new sample molecule velocity components after collision, the velocity 
components of the target molecule collision partner must first be found. Since there are 
two distributions of target molecules locally, f^ + and f^_, the distribution from which 
the target molecule collision partner is to be picked must first be decided. The fraction 
of all the collision partners that come from the f M+ distribution is given by X g A M+ - 
Then, for a picked random number, if R < X g A M+ > the collision partner velocity com- 
ponents are picked from the f M+ distribution. Otherwise, they are picked from the f M _ 
distribution. The equations for obtaining the velocity components of the target molecule 
collision partner are given in appendix C. 

If the V 2 t picked from the f^ + distribution is negative or, if the picked from 

f M- is positive, a collision does not occur because these target molecules are nonexis- 
tent, as discussed previously. In these cases, the sample molecule continues on from 
the point of no collision with its original velocity. 

After collision of the sample molecule with the chosen collision partner, the new 
components of velocity Vf i = 1, 2, 3, of the sample molecule (flow chart position 4) 
are given (ref. 7) by 



2 '/2 

f^ls + V lt> + V 1 - b > H 


(8a) 


V 2 s 4 <V 2 s + V 2t )+ i e<1 ' 2b2) 

* 1 9 !/2 

V 3s = \ <V 3s + V 3t> + V 1 " b > S 


(8b) 


(8c) 
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where is given by equation (B2) . The values of H and E are obtained by picking 
two random numbers which are used in the following equations: 


H = 2R r - f 

S = 2R„ - 1 l 


b 2 = H 2 + E 2 


( 9 ) 


o 

where b must be less than 1, or a new set of random numbers must be chosen to find 
H and E- 

The sample molecule history is then continued from the point of collision with its new 
velocity components. The sample history is completed when the sample molecule is re- 
turned to the surface from which it had been emitted. The numbers of sample molecule 
histories originating at the two walls were taken equal. 

The reason for emitting sample molecules from both walls was to improve the ac- 
curacy of the results. Especially for the cases involving small mean free paths, not 
enough of the sample molecules emitted from one wall reach the other wall to give suf- 
ficiently large sample sizes. This problem was avoided by emitting sample molecules 
from both walls. 


Scoring to Find Macroscopic Flow Properties 

The macroscopic fluid characteristics needed are density, temperature, velocity, 
and heat transfer across the channel. These properties are obtained at scoring positions 
located at various distances across the channel, as shown in figure 1 (flow chart posi- 
tion 5) . It some property of the molecule is designated by Q, the average quantity of Q 
transported per sample molecule across the scoring cross section p in the positive Xg 
direction can be written as 
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where S is the number of sample molecules passing across the scoring cross section 
P 

p in the positive Xg direction. Similarly, the average quantity Q transported per 
sample molecule in the negative Xg direction is 



( 11 ) 


Sp P_<v 2 > 


Since there is no net flow across the channel and all sample histories start and end 

at the same wall, S = S „ and 
’ +,P P 


(e + / v 2 f + d3y ) = -(<•- f v 2 f - d M 


( 12 ) 


Then, since p(V 2 Q) = P + ( V 2^) + P-( V 2^) > 


p + < v 2> S J 1 


(13) 


can be obtained. 

The number of sample molecules S + that pass the scoring cross section in the posi- 
tive Xg direction at p divided by the total number N of sample molecule histories 
leaving wall 0, can be related to the mass flux passing in the positive Xg direction at p 
by 


+,P 

N 


P + < V 2> 


P + < V 2> 


(14) 


wO 


where (p + (Vg) ) the mass flux leaving wall wO, is equal to P + w qC w o/ n 

V + /w0 


7 w0 

bining equations (13) and (14) results in 
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Thus, summing the inverse Vg velocity of the sample molecules as they cross the scor- 
ing positions yields p + and p_, the local densities of the molecule with positive and neg- 
ative components of Vg, respectively, and also the local density, p. These values are 
used in the assumed form of the target molecule velocity distribution in the next iteration. 
The average density across the channel is then obtained by averaging the density of all the 
scoring cross sections 


= J_ 

P + 0 Pf 



If Q is taken as Vj equation (15) becomes 


(p u i) 
P_ 

P+0 C w0 





(17) 


(18) 


Dividing by the appropriate densities obtained previously gives the local mean velocity 

for the molecules with Vg > 0, which is u^ + , the local mean velocity for the molecules 

with Vg < 0, which is u^_, and the local mean velocity u^. These values are also used 

in the target molecule velocity distribution in the next iteration. 

The kinetic energy per unit mass of a monatomic gas measured relative to the mean 

flow velocity is equal to the number of degrees of freedom times (1/2)R T. With this as 

& 
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the definition of temperature, the following equation can be written for three degrees of 
freedom: 


- R T = A (V’ 2 ) = A (V 2 > - A u 2 (19) 

2 g 2 2 2 1 


Then, from equation (15) 



The thermal velocity of the molecules with Vg > 0 is given by 



and similarly for the molecules with Vg < 0 



These values are used in the target molecule velocity distribution in the next iteration. 

The shear stress across the channel is given by P 12 = pCVjVg) . This can be shown 
to be equal to 


(p 12 ) p 


[p<VjV 2 >] 

p 


p + 0 C +0 

Ntt 1 / 2 




(23) 


The rate at which total energy is transferred across the channel is given by 
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(24) 


q = I[p(v(y 2 ))] 

T P 2 2 p 2 Nj7 1/2 


^ (V 2 ) - £ (V 2 ) 


o 

The transfer of the undirected kinetic energy across the channel is q = 1/2 p(V£(V' )} 
so that equation (24) becomes 


q * q T - u t P 12 


( 25 ) 


The macroscopic quantities of interest are found by using these equations. The 
values of p , p , C , C , u- , and u. are then compared to the assumed values used 
in the target molecule distribution, and the problem is iterated until they are in agree- 
ment. (The initial assumed values are taken from the collisionless solution. ) Also ob- 
tained at the same time from equations (16) to (25) are the local temperature, density, 
shear, and energy transfer. Twenty zones of equal size across the channel were used, 
and it was found that increasing this number had no effect on the results of the present 
calculations. 


RESULTS AND CONCLUSIONS 

First the case of stationary walls with a temperature ratio Vi of 0. 25 was treated. 
The local gas temperature divided by the temperature of the hot wall is shown in figure 3. 
The local density divided by the average density across the channel is given in figure 4. 
Figure 5 shows the heat flux transferred divided by the heat flux for the collisionless 
case, which is derived in appendix D. 

These results are in good agreement with the results of another Monte Carlo solution 
by Haviland (ref. 7). The linearized solution of Gross and Ziering (ref. 3) gives very 
good agreement with the present results except for the temperature profile. In the lin- 
earized solution of reference 3 the centerline temperature is taken midway between the 
wall temperatures. The present results have their centerline temperature significantly 
below this value for the Knudsen number of 2; hence the large disagreement in the results. 

The nonlinear results of Liu and Lees (ref. 4) of the temperature and density are of 
the correct magnitude but of a somewhat different gradient. The nonlinear heat-transfer 
results fall somewhat above the other solutions. 

Also shown are the slip continuum results. These results are obtained by using the 
fluid temperature near the wall (slip wall temperature) rather than the wall temperature 
in the continuum equations (appendix E). By using the slip wall temperatures that gave 
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Figure 3. - Temperature distributions. Wall-temperature ratio, 0.25; wall Mach number, 0. 
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Figure 4 - Density ratio for various Knudsen numbers. Wall-temperature ratio, 0.25; 
wall Mach number, 0. 
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Figure 5. - Dimensionless heat flux as function of Knudsen number. Wall-temperature ratio, 0. 25; 
wall Mach number, 0. 


the best fit to both the temperature and the density results of the present solution, the re- 
sults shown were obtained. However, when the heat transfer was calculated using these 
same slip wall temperatures and the continuum thermal conductivity the heat-transfer 
rate was considerably above that from the other solutions. Since the slip solutions with 
the continuum conductivity are only expected to apply close to continuum conditions, this 
result is to be expected. 

The heat-transfer rate at various positions across the channel is shown in figure 6. 
Theoretically this should be constant at each point across the channel (appendix E), al- 
though the present method does not impose this condition by its method of solution. This 
condition, however, is satisfied by the present results, as can be seen from this figure. 

To check the confidence limits of the present results, the solution was carried out 
for a Knudsen number of 0. 5 using 50 000 sample molecules. The temperature and den- 
sity had a value of 6, the 95 percent confidence interval, of close to 1 percent of their 
value, while the value of 6 for the heat-transfer rate across the channel was about 
2 percent of its value. The derivations and equation used in obtaining 6 are given in ap- 
pendix F. 

The second case treated is for that in which both walls are at the same temperature 
but one wall is moving at a wall Mach number, = u^ /C w q, of 4. The velocity pro- 
file is shown in figure 7. Since the profile is antisymmetric about x/D of 0. 5, the re- 
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Figure 6. - Dimensionless heat flux distribution for various Knudsen numbers. Wall- 
temperature ratio, 0.25; wall Mach number, 0. 
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Monte Carlo results, 
Knudsen number 
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Figure 8. - Temperature distribution for vari- 
ous Knudsen numbers. Wall Mach number, 
4; wall-temperature ratio, 1. 


suits are only shown for x/D from 0. 5 to 1. 0. The linearized solution (ref. 2) departs 
significantly from the present result obtained herein for the nonlinear large wall Mach 
number case. 

The temperature profile across the channel is shown in figure 8. There is a sub- 
stantial rise in temperature for larger Knudsen numbers because as the strongly directed 
stream of molecules leaves the moving wall it encounters the undirected stream. This 
gives a lower value for the average directed motion of the fluid. This decrease in di- 
rected kinetic energy is now included in the thermal motion giving the higher temperature. 

The shear stress for the channel is shown in figure 9. This has been nondimension- 
alized by the shear stress for the collisionless case (eq. (D5)). The results fall above 
the linearized solution of reference 1. The total energy transfer across the channel non- 
dim ensionalized by the total energy transfer for the collisionless case (eq. (D6)) gave the 
same values as the dimensionless shear stress for the cases studied. This is also shown 
in figure 10, where shear stress and total heat transfer at various positions across the 
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Monte Carlo results, 
Knudsen number 

□ lxlO 20 

O 2 

^ 1 


Closed symbols denote shear stress » , 

F p l, z/PaCwJ/^' 2 p l, 2/PA c wo) Kn -« 

Open symbols denote total energy transfer = 

p /2( lt/PA c wo)/( 7rl/2( lt/pA C S)) Kn , 00 



Nondimensional distance from wall, x/D 

Figure 10. - Dimensionless total energy transfer and shear stress. Wall-temperature 
ratio, 1; wall Mach number, 4. 
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channel are plotted. Again, from the laws of conservation and the moment equations, it 
can be shown that the shear stress and total heat transfer at various positions across the 
channel should be constant. Although this condition is not imposed by the present method 
of analysis, the results in figure 10 indicate that the shear stress and total energy trans- 
fer are nearly constant at the different positions across the channel. 


CONCLUDING REMARKS 

The Monte Carlo method gives reasonable results for the present problems without 
the necessity of extreme assumptions. The major drawback to the method is the increas- 
ing amount of computer time necessary to carry out the calculations as the Knudsen num- 
ber becomes smaller. 

Some typical running times for 50 000 sample molecules on an IBM 7094 are 40 min- 
utes for 1^ = 30, 50 minutes for K r = 2, and 75 minutes for K n = 0. 5. The large in- 
crease in running time as the Knudsen number is reduced shows that the procedure as 
presently formulated is not suited to very small Knudsen numbers. Generally, three it- 
erations of 50 000 samples were needed. 

The present method is very flexible and can readily be used in other transport prob- 
lems not readily solved by the more usual procedures. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, February 20, 1968, 

129-01-11-05-22. 
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APPENDIX A 


VELOCITIES OF MOLECULES LEAVING THE WALL 


The normalized Maxwellian distribution of the molecules at the wall wO with com- 
ponent Vg > 0 is given by 


*wo dv ' 3 = 


3 / 2^3 


"wO 


exp - 



dV 


,3 


(Al) 


(Note that integration of eq. (Al) over all velocity space with V2 > 0 will give unity. ) 
The distribution of velocities of the molecules leaving the wall wO per unit time per 

unit area is given by V 0 f /(V 0 ) . This distribution can be transformed to cylindri- 

Z w 1 +w0 

cal coordinates Vj = V r cos 6, V 2 = V 2 , Vg = V r sin 6 as 


V 2Vo dV ' 3 2V 2 V : 


<v 2 > 


+w0 


7TC. 


w 


exp 




"wO 


dO dV 2 dV r 


(A2) 


Equation (A2) can be written as the product of the following three single variables distri- 
bution functions: 



(A3a) 


(A 3b) 



(A3c) 


The velocity components of the sample molecule leaving the surface must be picked from 
these distributions. A convenient way of picking from a distribution for the high-speed 
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computer is to transform the distribution to a uniform distribution by setting the random 
number R equal to the cumulative distribution function. For example, 




f 0 , d0* = — 


27 T 


(A4) 


Then, by using a high-speed computer to generate a random number R between 0 
and 1, 6 can be obtained from equation (A 4) such that for a large number of samples 
picked in this manner the distribution in equation (A3a) will be satisfied. This technique 
is more fully discussed in reference 1. 

Similarly, V 2 can be picked from 



Or, since picking R is equivalent to picking 1 - R, 



(A6) 


so that by picking a random number 
ilarly 


, V 0 can be obtained from equation (A6) . 
V 2 ^ 


Sim- 




(A7) 


A similar procedure is used to generate sample molecules from wall wl using in- 
stead C wl and V 2 < 0. 
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APPENDIX B 


EVALUATION OF MEAN FREE PATH 


The incremental number of collisions per unit time d0 of a sample molecule moving 

3 

at velocity V through target molecules in velocity space increment dV f for hard 

S t 

sphere molecules is given, in reference 9, as 


P f t 3 

d 0 = — t y s dV7 

m R l 


(Bl) 


where V 0 is the magnitude of the velocity of the target molecules relative to the sample 

rv 

molecule velocity before collision, Vt, . = V. . - V_ .. The term V w can be written as 

Xv« 1 L. 1 S, 1 Xv 


V R = 


(V 


it - v i s ) 2 + ( y 2t - v 2«> 2 + - v 


1 1/2 


2s ; 


3t 


3s' 


(B2) 


If the local mean velocity u = Uj (which is along the Xj -direction) is subtracted from 
both the sample molecule and the target molecules, the collision rate for a Maxwellian 
distribution of target molecules can then be written as 




pS 


3/2„3 
mff ' C 


dV it dV 2t dV 3t 


(B3) 


In order to carry out the integration to obtain 0^, the reference coordinate system 
is rotated so that V* lies along the new 2-direction. This is accomplished by employ- 

S 

ing Eulerian angles (refs. 9 and 11) (see fig. 1(b)). The rotation can be considered to be 
a two-step process whereby the coordinate system is first rotated by the angle cp coun- 
terclockwise with axis along the 3-direction, where 


V 2s v is 

COS cp = — — , sm cp = - 


rs 


rs 


/ o 2 

and V rg = ( Vj g + Vgg ) . After this rotation of the 1, 2-reference plane, the plane 

containing the velocity vector V’ of the sample molecule and the new 2-direction is 

s 

perpendicular to the new 1-direction so that another rotation counterclockwise about an 
axis along the new 1 -direction by the angle 9 where 
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cos 9 = 


V 


T?’ S inO=-i?, + 


1/2 


places the 2 -direction along V’ . The matrices corresponding to the preceding trans- 

o 

formation and its inverse are 


where 


and 


where 




(B4a) 



^is 

V s 


V. 


2s 


V s 


3s 


V’ 

s 


■ V 3s V ls 

V rs V s 


• V 3s V 2s 

V rs V s 



(B4b) 


^ = A ij V it 


(B5a) 



■ V is 


rs 


'2s 


V s 


- V 2s V 3s 

V rs V s 



(B5b) 


Then equation (B3) becomes 
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rmr 3 / 2 C 3 


exp 


V t ,2 \r 2 2 21 1 / 2 s 

“ J fit + < V 2t " V s> + V 3t J dV t’ 


(B6) 


Transforming into spherical coordinates 


V it = V t cos 9 sin <f>t V 2t = V t cos <p> V 3t = V t sin 9 sin ^ 


(B7a) 


results in 


= [exp(- M 3 )] /i t 2 M R sin cp dcp dd dp t 

mn' i ' 


(B7b) 


where /i is the nondimensionalized velocity V'/C and where 


^R = (^t + ^s “ 2 ^s cos <P) 


\ 1/2 


(B8) 


To obtain the total collision rate for the target molecules over all velocities, equa- 
tion (B7b) must be integrated over cp, 6, and Integrating over 0 from 0 to 2ir and 
over cp from 0 to it gives 


[sxp(-m 2 )Jm 2 M r> a d ^t 

mn ' 


/-7T 

where M R ^ ~ Jq M R (sin cp/2)dcp 


or 


M 2 

^s+rf * ^ s >^t 


m R, A "A 




Mt+ £ * 


(B9) 


(BIO) 


Integrating equation (B9) over from 0 to °° yields the total collision frequency 
© M for a sample molecule moving at velocity ju g through a Maxwellian gas: 
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(Bll) 


7 M 


pSC 

exp 

m 

r-i 

1 


+ (erf fx g ) Mg 



The nondimensional mean free path is then found, as discussed in reference 10, by 
dividing the sample molecule velocity by the total collision rate and using the hard- sphere 
definition of Knudsen number Kn to give for the dimensionless mean free path 



If there were two different Maxwellian gases in the incremental volume in space, 

distributed over velocity space by the functions f Ma and f M ^, respectively, the total 

collision rate would be given by the collision rate for the molecules from f Ma plus the 

collision rate for the molecules from f M ^, 0 = © Mft + © M £, where © MQf and © M ^ are 

given as in equation (Bll) with appropriate densities and mean thermal velocities. Since 

X o is given by V /©, 

& s 


x 


s 


1 _ 1 
0 Ma t e M/3 

V c V c X M« X MP 

o b 


(B13) 
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APPENDIX C 


TARGET-MOLECULE COLLISION PARTNER 


The distribution of target-molecule collision partners in velocity space for a sample 
molecule passing through target molecules in a Maxwellian distribution is given by divid- 
ing equation (B7b) by equation (Bll): 


d© 


M 


exp^Mt 2 ) 


p^ Mr sin cp d cp d& dp t 


0 


M 


Jr[exp(-M 2 )] + 7r 3 / 2 (erf p s )hz g +~- 


(Cl) 


The distribution in Q for the target-molecule collision partners is readily seen to be 


f dO = M 
6 2n 


(C2) 


where the value of 6 can be picked from this distribution as discussed in appendix A by 
employing the relation 


e = 2?7R 


6 


(C3) 


The distribution of p t for the target molecule collision partners is obtained from 
the marginal distribution 


f dp 


t 


4 [ eXP (-P t )lp t pR, A dM t 


exp (-p 2 ) +>r 


1/2 (erf P s )(p g + 


s/ 


(C4) 


The value of p t from this distribution is found from 

y(p s ,p t ) + e (P s >P t ) 


%" 


exp(-p 2 )] + >r 1/2 (erf P s )(p s + -- 


2/x, 


(C5) 
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where 


3 

y(M s > M t ) = -2M s M t [exp(-M t 2 )] + M s (erf M t ) - | [exp(-^J)J 




H 


[exp^t 2 )] 


2 ^ +— (erf (i t ) 

Ms 1 - ' " J 2 M s 


>when M s > M^ (C6) 


e = 0 




and 


y(Mg, M t ) = y(M s , m s ) 


- 2 [exp(-M^|fe + i + M? 


, 2N v<m 2 V whe " , ‘ s< ' lt <C7) 

2 exp (' M s)J( - r + 1 ' 


To find (p for the target molecule, the distribution of target molecules can be writ- 
ten as a product of a marginal times a conditional distribution 


d°(M t ,<P) 

= f(ju t )f(«p|u t ) 


© 


(C8) 


Then cp must be found from the conditional distribution 


f ifp | M t ) = 


sin (p d <p 


2ju 


R, A 


(C9) 


The value of <p for a given is found from 


R„ = 
<P 


4- 4 - ^ 


/ 2 2 \ 3/2 

(p t + M s + 2 M t M s ; -|M t -M 


(CIO) 
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where 


4=^ 


o 

+ M s " 2 M t M s cos 


\3/2 


After cp, and 9 are found, the velocity components V^, V^., and 
by using the relations given by equation (B7a). Then, from equation (B4) 


V 


It 




V 3s V is\ 

, V rs V s/ 


+ u 


1 


V 


2t 



- V 3t 


V 2s V 3s\ 

V V’ / 
As s/ 


V 


3t 




are found 


(Clla) 


(Cl lb) 


(Cllc) 
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APPENDIX D 


COLLISIONLESS SOLUTION 


For the case of negligible molecular collisions the local distribution in the channel 
is given by 


P f oo 



v: 


,2 


exp - 


wO 


exp - 


"'wO/ 




"wl/ 




V 2 >0 




V 2 <0 


(Dl) 


where p + Q is the density of molecules with V 2 > 0, while p ^ is the density of mole- 
cules with V 2 < 0. The local density is then given by 


P(V 2 Q) / =P( 1) 

q-i/v 2 


= P=P + 0 + P-l 

Since there is no net flow across the channel, 

2> =-^(P + 0 C w0-P-l C W l)= 0 

It ' 

From the first relation 



where, from the second relation 


(D2a) 
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1+ 5w» i+t; 1 / 2 


The flow parallel to the channel wall is given by 


U 1 = (Vj) =-±- Ul mn+ ^i Ul 


7 1-wO T 71.wl 


The local temperature as can be seen from equation (19) is given by 


_ T _ = _ JL _ (y 2 ) - - M 2 
T w0 3 C 2 q 3 


which when evaluated gives 


t = ? /_±2 M 2 + M 2 , - M 
3 \ p w0 p wl 




The shear stress (eq. (23)) P 12 = P( v i V 2) can be found to be 


77 P 12 p + 

~^r~ - -f < M wo - “wi* 
p2R g T wO P 

The total energy transferred is given by equation (25) as q T = p/2 (V ,2 V 2 ) + UjPjg. 
For the collisionless case, the relation becomes 


1/2 2 2 
77 q T _ p + 0 _ p - 1 t 3/2 + P+0 M w0 _ M M wl t l/2 

p(2R g T wO )3/2 P P ^ P 2 P 2 wl 
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APPENDIX E 


CONTINUUM SOLUTIONS 

Maxwell’s equation of transfer for the steady- state condition can be written 

— (p(V.Q)) = m A(Q) (El) 

dx i 1 

where the convention is used that repeated subscripts are to be summed over. Here, 

A(Q) denotes the rate of increase in some molecular quantity Q due to molecular col- 
lisions. This term is zero for cases where Q is conserved during collision. This is 
true when Q is equal to 1, V, or V^, corresponding to conservation of mass, momen- 
tum, or energy, respectively: 

Q= 1: 


9X i 


(puj) = 0 


(E2) 


Q = V: 


o = -J-(p<v V.» 

3Y 1 j 


=^[p< v i +u i ,<v j +u J»] 


■=— <P<V!V!> + puu) 
PIY 1 J 1 J 


9X i 


(p(VjVj) ) + pUj 


3u. 


+ u. 


9X i 


3(pUj) 


9X i 


Inasmuch as p(V|V p = Py and 9(pu i )/9X i = 0 by equation (E2), the preceding equation 
can be written as 
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(E3) 


9P.. 9u. 

— ii + pm — 1 = 0 

aXj 1 ax t 


Q = V 2 : 



9X i 


[p<(V! + u.)(V' + u) 2 >] 


= -1- (p(V|V' 2 >) + -L (u.P. i ) + -1- (2u i P ii ) 
3Xj ' 1 ' 3X. 1 " 9X- 1 1J 


9x i 



By equations (E2) and (E3) the sum of the last two terms of the preceding gives 
(2 P a 9u j/ 9X.), so that 



+ — (UjPjj) + 2 Pii 
9X t 1 1] 


9u. 


ax. 


= o 


(E4) 


In addition, there exist the two well-known relations for a fluid, namely, 

q s 1 p(V!V' 2 ) = ~/c — 

1 2 1 ax. 


(E5) 


P ii 


= 6 ijP 


2p 


v 9X i 



2 dX i 




(E6) 


Equations (E5) and (E6) serve to introduce the gas conductivity k , the viscosity p, 
and the static pressure P into the relations. In this context, note that the introduction 
of equation (E6) into equation (E4) yields the Navier- Stokes equations. 

With respect to the problem treated herein, where the only changes occur in the Xg- 
or across-channel-direction, equations (E2) to (E4) become 
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(E2): 


3u 


3X 


1 = 0 or u 1 = g(x 2 ) 


(E7a) 


(E3): 


3P 2i 3 Uh 

i =0 or P 1<? = -ju = Constant 

3X 2 3X 2 


(E7b) 


where equation (E6) has also been employed, 
(E4) : 


dXr 


(q 2 + P 2lUl ) = 0 


so that by equation (2 5) 


q T = q 2 + MjPj 2 = Constant 


(E7c) 


For hard- sphere molecules, the temperature- dependent conductivity and viscosity 
can be written for the continuum case as (ref. 6) 


K = K Q t 


1/2 


where 


K 0 = 


75m K g^w0 


64a 2 (27 r) 1 / 2 


(E8) 




and 


ju. = /u n t 


1/2 
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where 


5mCw n 

Mo= — i — I7i <E9) 

16(7 (2jt) A/ ^ 

For the ease where both walls are stationary, equation (E7c) becomes 

q 9 = = Constant (E10) 

* D 8X2 

where the dimensionless coordinate Xg = Xg /D has been introduced. Substituting equa- 
tion (E8) into equation (10) and integrating give 


1=0 ■ [ l ' (j 


} - 


Because the pressure across the channel is uniform, the ideal gas law becomes 

pt M=0 = Constant 

Then, employing the relation p^ = J p dXg shows that 


1 + + fc wl 


3t M=0 


Also the heat-transfer rate across the channel is 


fc 3T K 0 t T w0 at M=0 


M=0 D 3Xci 


which upon evaluation becomes 


p A C wO, 




M=0 
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For the case where wall wl is moving at velocity ju equation (25) applies: 

q T = q 2 + u iP 12 


_ K 3T _ n du l „ 
D 0Xg D SXg 


K 


_0T + 1 _^1 

CpD y p 9Xg k 2 9x 2> 


Inasmuch as c = 5/2 R and Pr = nc /k = 2/3 for a gas, this equation becomes 

lr O ir 


[c T 


0x 2 \ P 3 


u i \ c„ 

+ — |= -^Dq n 


The right side can also be written as a partial derivative with respect to x 2 by 
making use in turn of the relations M c p A = 2/3 and equations (E7b) and (E7c): 


' _ED,T= ‘(j i ,D “ T 


2 q T 


3u 1 


3 P 12 3x 2 


d 1 2 q T 


u, 


9x 2 \ 3 P 12 


Thus, equation (E14) is transformed to 


3 fcT + !i- 3 ^u.l=0 
Sx 2 \ P 33 P 12 


Dividing through by c p T w0 = 5/2 R g T w0 = 5/4 CT Q yields 


(E14) 
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8 


3x r 


M = 0 


15 


15 P 12 C wO 


For the case where both walls are at the same temperature, that is, t w g = t^ 
integration of equation (El 5) gives 


t = -2- M(M - - M) + 1 
15 W1 


and the relation 


q T M wl P 12 C wO 


If equation (E9) is substituted into equation (E7b), the result is 


12 


u t V 2 c 
M 0 L L w0 _3M 

D dx 0 


= Constant 


Substitution of equation (E16) yields 



- M) + 1 = Constant 

0x 2 


Integrating the preceding equation gives 


x 2 = 



- M wi + ( 15) 1 / 2 



r 



. -i 

sin 


M , - 2M 
wl 

(l5 + 



+ 15<> 


sin - * 


M 


wl 




1/2 



(El 5) 

= 1, 

(E16) 
(E 16a) 

(E17) 


(E18) 
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and the relation 


r 


, 1 / 2 , 


12 5 n 


p A C wO 


32 


Kn 1 M wl 





(E19) 


Equation (E18) gives the velocity profile across the channel and may be used in conjunc- 
tion with equation (E16) to give the temperature distribution across the channel. The 
shear stress is given by equation (E19) and the total heat flow by equation (E16a). 
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APPENDIX F 


CONFIDENCE LIMITS OF RESULTS 


Following the derivation in reference (1), the confidence limits on the present re- 
sults can be obtained as follows. For some number of sample molecules a quantity g is 
obtained. Repeating this process k times will give a sample set of g values. The av- 
erage of the k values would give g, which is the average value of the sample set. If an 
infinite number of values of g were obtained, the average value of g would be the true 
average denoted as (g) . It is necessary to determine a value of 6 such that the absolute 
value of the difference between (g) and g is less than 6 with a 95 percent probability. 
The central limit theorem states that the distribution of sample means f(g) approaches a 
normal distribution with a mean value equal to the population mean (g) and with a vari- 
ance equal to the population variance divided by k, the number of measurements used to 
determine g. Since the 95 percent confidence interval 6 corresponds to twice the stan- 
dard deviation for the average of k values, 6 = 2S/(k) ' , where S is the approxima- 
tion to the variance of the population obtained by finding the variance of the sample set 


,2 _ [g 2 - fe> 2 ] k 

k - 1 


(FI) 


This gives the 95 percent confidence interval 6 as 


[<g> 



(F2) 
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